Calculations of canonical averages from the grand canonical ensemble 
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Grand canonical and canonical ensembles become equivalent in the thermodynamic limit, but 
when the system size is finite the results obtained in the two ensembles deviate from each other. In 
many important cases, the canonical ensemble provides an appropriate physical description but it 
is often much easier to perform the calculations in the corresponding grand canonical ensemble. We 
present a method to compute averages in canonical ensemble based on calculations of the expecta- 
tion values in grand canonical ensemble. The number of particles, which is fixed in the canonical 
ensemble, is not necessarily the same as the average number of particles in the grand canonical 
ensemble. 
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I. INTRODUCTION 



Whenever we need to work with a system of many 
quantum particles it is much easier to perform the cal- 
culation in the grand canonical ensemble than the cor- 
responding calculations in the canonical ensemble. For 
example, the calculations of the grand canonical parti- 
tion function for ideal gas of fermions or bosons is triv- 
ial, whereas the computation of the canonical partition 
function for the the same system becomes a formidable 
task even for a small number of particles. Provided that 
the system is thermodynamically large, the canonical 
and grand canonical descriptions agree with each other. 
There are many quantum systems where the canonical 
description is more appropriate, these include hot nuclei 
PI, ultrasmall metallic grains Q, Bose and Fermi gases in 
atomic traps @,IU and atoms in plasmas Therefore, it 
is important to have a practical theoretical and computa- 
tional method which enables us to extract canonical av- 
erages from corresponding grand canonical calculations. 

The problem, which we would like to solve in this pa- 
per is the following. Suppose that we can perform calcu- 
lations (or measurements) in grand canonical ensemble 
which is characterized by temperature T and chemical 
potential fi. We would like to compute expectation val- 
ues of an arbitrary quantity O in the canonical ensemble 
with temperature T and the number of particles n by 
using only the averages in the grand canonical ensemble. 
The number of particles n, which is fixed in canonical en- 
semble, is not necessarily the same as the average number 
of particles (N) in the grand canonical ensemble. 

Earlier works on particle number projection in grand 
canonical ensembles have been performed to treat hot 
nuclei [H, 10,0, Bose-Einstein condensation! 3L Uj, Kj| and 
to formulate canonical statistical mean field approxima- 
tion for mesoscopic systems jlOj. Our approach is very 
different. First, we do not rely on projection operators 
but extract information from the grand canonical aver- 
ages by inverting fluctuation matrix. Second, (N) does 
not necessarily equal to n although it can be. 



II. THEORY 

A quantum mechanical system has the Hamiltonian 
H. The Hamiltonian H commutes with the number of 
particles operator N 



[H,N} = 0. 



(1) 



Therefore the Hamiltonian and the number of particles 
operator have the same set of eigenvectors: 



N\na) = n\na), 



H\na) — E na \na). 



(2) 



(3) 



Let O be an arbitrary operator. The average in the grand 
canonical ensemble is 

(0) = ^Y, e ~ PiE,ia ^ n) ( na \ \ na )> ( 4 ) 



with Z being the grand canonical partition function 

z = Y^e-^ B ™~» n \ (5) 



and [3 = 1/ kT. The average in canonical ensemble is 

(0>„ = ^5>-^<HO|na>, (6) 

a 

with Z n being the canonical partition function 

Z n = J2e-^». (7) 

a 

We can consider (0) n as a function of n and expand it 
as the power series: 



(0) n =^q k n k . 



(8) 



fc=0 
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With this expansion the grand canonical average (|4|) be- 
comes (appendix A): 



(0) = C£q k N k ). 



(9) 



fc=0 



Next, we take canonical expectation value and 
add/subtract the grand canonical average 



(0) n = J^q k n k + (0)-(^q k N k ). 



(10) 



k=0 



k=0 



Then we cancel k — terms in the both sums and the 
canonical expectation value becomes 



(0) n = (0)-J2q k ((N k )-n k ). 



(11) 



fc=i 



This expression for canonical average of the operator O 
involves only grand canonical expectation values. The 
coefficients q k are yet to be calculated. To determine 
them we introduce the operator 



= 0-J2q k N k 



(12) 



fe=i 



and compute the expectation value (Of(N)), where /(A) 
is an arbitrary function of A. We show in appendix A 
that this expectation value can be split: 



(Of(N)) = (0)(f(N)). 



(13) 



Since /(A) is an arbitrary function of N, Eq.fT3|) is 
equivalent to the following system of equations 



(ON) = (0)(N) 
(ON 2 ) = (0)(N 2 ) 



(ON k ) = (0)(N k ) 



(14) 



If we use the explicit form for operator O (fT2")) , the syste 
of equations (fT4"]) becomes: 



an 



J2 IkAkm = (ON r 



(0)(N T 



(15) 



fc=i 



where the matrix A km is built from the fluctuations 

A km = (N k+m ) - (N k )(N m ). (16) 

The system of linear algebraic equations (|15p along with 
the expansion pip represent the main result of the paper. 
The similar equations (although for the case (N) = n) 
were obtained by the methods of thermofield dynam- 
ics in the context of superconducting nuclei at finite 
temperature[i"ll li~2| . 



It can be readily shown by the direct differentiation of 
the partition function ([5]) that 



1 



Z 



gk+m Z d k Zd m Z 



(17) 



Now we prove the convergence of the expansion (jlljl . 
It is possible to demonstrate based on simple consider- 
ation that the difference between canonical and grand 
canonical averages is [l3[ 



(O) (0) r 



1 



(A)' 



(18) 



The calculation in appendix C shows that for k > 1 
q k ~l/(N) k . 



(19) 



Suppose that n = (N). Then {(N k ) - n k ) - (N) k ^ 1 
for large (N). Therefore, the corrections to the grand 
canonical average in Eq. fTTj) become in this case 



oo . oo 



(20) 



k=l 



k=l 



where c k are some coefficients. Comparing (fT5| and (|20|) 
we see that ^2^ =1 c k is finite. It means that the expan- 
sion (fTTj) converges when (N) = n. We would like to note 
at this point that each term in the sum pip is ~ 1/ (A) . 
It means that the method becomes computationally ef- 
ficient when it is applied to the large systems, since one 
needs to include less terms in the expansion to achieve 
the same accuracy in this case. Let us consider the cor- 
rections to the grand canonical average in Eq. P"Tj) for the 
case n = (A) + j where j is some integral number 



(0)-(0) n = Y,1k((N k )-[(N)+j} k ). 



(21) 



k=l 



We assume that j/(N) < 1. If we substitute the ex- 
pression for (A fc ) (|C17p into (|2"Tj) and use the binomial 
expansion up to the first order in j / (N) for [(A) + j] k 7 
we obtain 



(0)-(0) n = J2Qk(N) 
fc=i 



fc _ 1 cfc(/t- 1) 



1 - 



■2j 



c(k - 1) 



(22) 

The series (|2"2"| always converges for j = as we just 
demonstrated. To prove the convergence for j ^ we 
split the sum (j2"2")) into two parts: the first part is for 
< k < K and the second part is for K < k < oo, 
where K is some positive integer. The first part is al- 
ways finite and ~ 1/(A). By choosing K, we can always 



make 



1 



2j 



c(fc-l) 



1 for k > K. Therefore, the conver- 
gence of the expansion (fTTj) for (A) = n case implies the 
convergence for any finite difference |(A) — n\ for which 
^ N (n)^ ^ 1- W e note that these arguments may not 
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TABLE I: Occupation numbers and total energy for fermions. 
FD refers to the Fermi-Dirac statistics in grand canonical en- 
semble ({N) = 4). The number of particles in the canonical 
ensemble is 2. k max is the number of terms included in ex- 
pansion (lllfl . 













kmax 






1 


Si 


FD 


1 


2 


3 


4 6 


exact 


1 


1.0 


0.98 


0.92 


0.87 


0.86 


0.87 0.87 


0.87 


2 


2.0 


0.95 


0.80 


0.69 


0.67 


0.67 0.68 


0.68 


3 


3.0 


0.87 


0.52 


0.34 


0.34 


0.30 0.30 


0.30 


4 


4.0 


0.72 


0.07 


0.01 


0.06 


0.11 0.12 


0.12 


5 


5.0 


0.48 


-0.31 


0.10 


0.06 


0.04 0.04 


0.04 


total 


energy 


10.77 


2.82 


3.77 


3.78 


3.79 3.79 


3.79 



TABLE II: Occupation numbers and total energy for bosons. 
BE refers to the Bose-Einstein statistics in grand canonical 
ensemble ({N) = 4). The number of particles in the canon- 
ical ensemble is 2. kmax is the number of terms included in 
expansion {TTJ. 













kmax 






1 


£l 


BE 


1 


3 


5 


7 11 


exact 


1 


1.0 


3.43 


1.52 


1.52 


1.47 


1.44 1.42 


1.42 


2 


2.0 


0.40 


0.33 


0.33 


0.37 0.39 0.39 


0.39 


3 


3.0 


0.12 


0.10 


0.10 


0.11 


0.12 0.13 


0.13 


4 


4.0 


0.04 


0.04 


0.04 


0.04 


0.04 0.05 


0.05 


5 


5.0 


0.01 


0.01 


0.01 


0.01 


0.02 0.02 


0.02 


total 


energy 


4.81 


2.68 


2.69 


2.77 


2.82 2.84 


2.85 



work when the matrix A nm is singular. For example, in 
the case of low temperature Fermi gas, all matrix ele- 
ments A nrn tend to zero, therefore canonical and grand 
canonical descriptions may deviate from each other in the 
thermodynamic limit due to the persistent existence of a 
few-particle fluctuations in the grand canonical ensemble 

III. EXAMPLE CALCULATIONS 

As an example we consider the system of noninteract- 
ing quantum particles distributed on ni eve i s single parti- 
cle energy levels with energies ei . The logarithm of the 
grand canonical partition function is 

n levels 

\nZ = i 1 ]T ln(l + 7]e fJ ^- £,) y (23) 
1=1 

where r\ = +1 is for fermions and -q = —1 is for bosons. 
We set = 1, ei = I, and ni eve i s = 5 in all our calcula- 
tions. 

First, we extract averages in canonical ensemble for 
the smaller system from grand canonical ensemble for 
the larger system. We select the chemical potential ^ in 
such a way that the average number of particles (N) in 



TABLE III: Occupation numbers and total energy for 
fermions. FD refers to the Fermi-Diract statistics in grand 
canonical ensemble ((N) = 2). The number of particles in 
the canonical ensemble is 4. kmax is the number of terms 
included in expansion (|11[) . 



1 


£l 


FD 1 


2 


3 


4 


5 


exact 


1 


1.0 


0.80 1.18 


0.92 


0.98 


0.99 


0.99 


0.99 


2 


2.0 


0.60 1.18 


1.01 


0.97 


0.93 


0.97 


0.97 


3 


3.0 


0.36 0.91 


1.02 


0.96 


0.97 


0.91 


0.91 


4 


4.0 


0.17 0.51 


0.70 


0.72 


0.73 


0.77 


0.77 


5 


5.0 


0.07 0.23 


0.35 


0.38 


0.37 


0.36 


0.36 


total 


energy 


4.10 9.42 


10.53 


10.54 


10.55 


10.55 


10.55 



the grand canonical ensemble is 4. We would like to ex- 
tract the information about the canonical ensemble with 
n — 2 particles from this grand canonical ensemble. We 
compute the occupation numbers and then all physical 
quantities like total energy can be calculated with the use 
of these occupation numbers. To start our calculations 
we set O = ni, where ni is the operator of the number 
of particles on level I. Then we solve the linear system 
of linear equations (fl~5|) to find the coefficients qk and 
use these qkS to calculate the grand canonical occupation 
numbers by Eq. pip . The matrix elements A nm are com- 
puted by Eq. (|17|) and with the help of recurrent relation 
from appendix B. The matrix element in the right hand 
side of Eq. (fT5|) is computed as the following derivative of 
the grand canonical partition function (|23| : 

The results of these calculations are shown in Table 1 
(fermions) and Table II (bosons). The fermionic and 
bosonic systems both show the convergence to the exact 
results as we include more terms in expansion (fTTj) . The 
convergence for bosons is slower than that for fermions. 
It is due to the fact that the fluctuations of the occupa- 
tion numbers (An^) = (ni) — r)(ni) 2 tend to be larger for 
bosons (rj — —1) than for fermions (77 = +1). 

The method also works in the opposite direction, there- 
fore we can compute averages in canonical ensemble for 
the larger system using grand canonical averages for the 
smaller system. We select the grand canonical ensemble 
with (N) = 2 and we compute the occupation numbers 
in the canonical ensemble of n = 4 particles. Table III 
shows the results of these calculations for noninteracting 
fermions. The convergence to the exact values is as good 
as in the previous case, thereby it demonstrates that the 
method can be also used to extract the canonical en- 
semble information for the larger system from the grand 
canonical calculations of the smaller system. The very 
similar results were obtained for bosons and we do not 
show it here. 
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IV. CONCLUSIONS 
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APPENDIX A: USEFUL IDENTITIES 



I 

Proof that if (0)„ - EZo1kn k then (O) = (EZ 1kN k ). 

(O) = ^^ iEna ^ n) (na\0\na) = |E^"E 

na n a 

1 oo 1 OO 

= I £ e^Z n £] q k n k = \ £ e -«*»-/«0 £ ftn * = I £ 



k=Q 



0\na) = i^^(0), 

n 

oo oo 

J2<lkN k \na) = (Y / qkN k 



k=0 



k=0 



k=0 



Proof that (Of(N)} = (0)(f{N)). 



(Of(N)) = i^e-^-^(HO/(iV)|na) = ± E e-^ E ^ n \na\0\na)f(n) 

na na 

-Ee^ n /(n)^e-^(na|0|na) = - E e^ n f(n) E e"^- (na|0 - E %iV fe |na} 



| E e^V(n)^« {{0) n - £ 5fe n fc ) = i £ 

n \ fe=l / n 



fc=l 



e^ n f(n)Z n q = (f(N))q Q = (f(N))(0) 



APPENDIX B: RECURRENT RELATION FOR 
THE CALCULATIONS OF THE DERIVATIVES 

We define d m = d m /(d^) m . 
Let 



Then 



Therefore 



* = InZ, Z m = d m Z. 



Z m = d m {e v ). 



Z m+1 = <9 m <9(e*) = <9 m (*'e*) 



E 

k=0 



k\{m — k)\ 



z ^(m-fe+l) 



where 



^r(fc) _ gkjg 



(Bl) 



(B2) 



(B3) 



(B4) 



Explicitly, 

^(°) = * = 7 ? E lnB ^ S / = (l + 77exp{/3(/i- £/ )})- 

(B5) 



r/ = 1 for fermions and = — 1 for bosons. 
Assume that 

* W = "E E^o- 

(=1 cr=0 

Since 



(B6) 



9(50-" = -vP{(Bi)- a - (Bi)-"- 1 }, (B7) 
we get then 

^ leve ls k 

= = E E -'^{(^rM^i) - " -1 }- 

(B8) 
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Since, according to (|B6 



Thus 



*<*+i) = £ £ a fe+i (B0 - (B9) 

;=i er=0 



we find that 



Evidently, 



-/3(<ra* - (a - l)a*_ x ). (BIO) 



oJ =j3, a\ = -p. (Bll) 



APPENDIX C: SCALING WITH THE NUMBER 
OF PARTICLES 

We use notations (B1)-(B5) from appendix B. In these 
notations 



(N 



m\ _ Z m 



Z0" 



(CI) 



Plugging this into Eq. (|B3p we get: 

m 

(N m+1 ) = E 



fc=0 



k\(m — k)\ 



(N k ) nm _ k+1 ■ (C2) 



In the thermodynamic limit, the sum over I in (|B5[) can be 
replaced by the integral. Since the corresponding density 
of states is proportional to the system volume V , e.g. in 
three dimensional space it becomes (fi is the mass of the 
particle) 



(C3) 



we see that - 7. Since (iV m ) must be proportional 
to V m , the term with k — m in Eq. (|C2[) gives the leading 
contribution. Retaining in (|C2[) the two leading terms, 
we get: 

*(2) 

(jV m+1 ) « (iV m ) (AT) + m— (iV™- 1 ) (C4) 



with 



(N*)-{N) 2 = J2((ni)-r){nif)- (C5) 



Here (nj) are the BE or FD occupation numbers. Since 
(A 7 ) > 1, we assume in Eq.([C4f that 



T 2 



(N m ) = (N) m + ^—a m (N) m - 2 + 0((N) m -' 1 ). (C6) 



Inserting (1C6|) into (|C4j) and retaining the leading terms, 
we get 



a m + m, a m = m(m — l)/2. (C7) 



(iV m ) « (iV) m + 2^ W(m , 1} <7V) m - 2 . (C8) 
p- 2 2 

Now, we consider the terms (niN m ). Following Eq. (f2"4]) 
we differentiate Eq. (|B3j) over ei, and retaining only the 
leading terms in V, we get: 

{ ni N m+1 ) = ( ni N m ) (N) + (N m ) ({m) - r) { ni ) 2 ). (C9) 

We shall look for the solution in the form 

( ni N m+1 ) = ( ni ) (N) m+1 



+7m+i (Nr+oaNr- 1 ). 



(CIO) 



Inserting this equation into (|C9[) . we obtain 7 m +i = 7-, 
(ni) - rj (m) and 

{ ni N m+1 ) a (n,) (7V) m+1 



-(m+l)((n/)-77(7i ; ) 2 )(iV)' 



(Cll) 



Next, we consider the system of linear equations (|15p . In 
the thermodynamic limit, we retain only the first terms 
in Eqs. ljClip and (|C8p . If we apply a "weak" thermody- 
namic limit and retain the next-order terms in Eqs. (|Cll|) 
and |C8l) . then 



(C12) 



(n,iV) r 



( ni ) (N) m a m (( nj ) - V (nif) (N)™- 1 (C13) 



Therefore, assuming that *( m ) ~ V ~ (AT) we get 



(AT)' 



ft (JV> 



fc+ro— 1 



(NY 



therefore 



q k ~l/(N) h 



(C14) 



(C15) 



(C16) 



The same (N) dependence can be also obtained if 
we express the particle number operator ./V in terms 
of creation/annihilation operators and apply Wick's 
theorempjl to matrix elements (7V m ) and (niN m ). 

Using the fact that *( 2 ) ~ (N) we transform Eq. (|C8)) 
to the form 



(N m ) a (N) m + 



(C17) 



Here c is some constant which does not depend on m and 
(AO- 

It has not escaped our notice that Eqs. (|C4[) and (|C9[) 
break down for bosons at critical and lower temperatures, 
since, due to the Bose condensation, ty( m > becomes pro- 
portional to V m . 



6 



[1] R. Rossignoli and P. Ring, Ann. Phys. 235, 350 (1994). 
[2] J. von Delft and D. C. Ralph, Phys. Rep. 345, 61 (2001). 
[3] H. D. Politzer, Phys. Rev. A 54, 5048 (1996). 
[4] C. Herzog and M. Olshanh, Phys. Rev. A 55, 3254 (1997). 
[5] F. Gilleron and J.-C. Pain, Phys. Rev. E 69, 056117 

(2004) . 

[6] Tomas R. Rodriguez, J. L. Egido, and L. M. Robledo, 

Phys. Rev. C 72, 064303 (2005). 
[7] K. Tanabe and H. Nakada, Phys. Rev. C 71, 024314 

(2005) . 

[8] H. Nakada and K. Tanabe, Phys. Rev. C 74, 061301(R) 

(2006) . 

[9] I. Fujiwara, D. ter Haar and H. Wergeland, J. Stat. Phys. 



2, 329 (1970). 

[10] S. A. Ponomarenko, M. E. Sherrill, D. P. Kilcrease and 
G. Csanak, J. Phys. A: Math. Gen. 39, L499 (2006). 

[11] D. S. Kosov and A. I. Vdovin, Z. Phys. A 355, 17 (1996). 

[12] D. S. Kosov and A. I. Vdovin, Izv. RAN (ser. fiz.) 60, 94 
(1996). 

[13] D. A McQuarrie, Statistical Mechanics, University Sci- 
ence Books, (2000). 

[14] S. P. Bowen, Y. Zhou, and J. D. Mancini, Phys. Rev. B 
46, 1338 (1992). 

[15] G. C. Wick, Phys. Rev. 80, 268 (1950). 



